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Cell-to-cell variability in cell death: can systems 
biology help us make sense of it all? 

X Xia 1 ' 2 ' 3 , MS Owen 1 ' 2 ' 3 , REC Lee 1 ' 2 and S Gaudet*' 1 ' 2 

One of the most common observations in cell death assays is that not all cells die at the same time, or at the same treatment 
dose. Here, using the perspective of the systems biology of apoptosis and the context of cancer treatment, we discuss possible 
sources of this cell-to-cell variability as well as its implications for quantitative measurements and computational models of cell 
death. Many different factors, both within and outside of the apoptosis signaling networks, have been correlated with the variable 
responses to various death-inducing treatments. Systems biology models offer us the opportunity to take a more synoptic view 
of the cell death process to identify multifactorial determinants of the cell death decision. Finally, with an eye toward 'systems 
pharmacology', we discuss how leveraging this new understanding should help us develop combination treatment strategies to 
compel cancer cells toward apoptosis by manipulating either the biochemical state of cancer cells or the dynamics of signal 
transduction. 
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Facts 

• Responses to apoptosis-inducing drugs or ligands are 
often heterogeneous. 

• Single-cell measurements are key to accurate character- 
ization and modeling of heterogeneous responses. 

• State-space maps allow contextualization of protein func- 
tion in cellular response. 

• Single-cell signaling dynamics can encode information on 
cellular responses. 



Questions 

• Can computational models of cell death incorporate multi- 
ple sources of cell-to-cell variability? 

• Can multivariate analysis approaches direct the design of 
effective combination therapies? 

• Can we improve on current anticancer therapies by 
optimizing treatment sequence and dosing schedules? 

Cell death in its various forms, including apoptosis and 
necrosis, has an essential role in development and in tissue 
homeostasis. When cell death regulation is derailed, disease 
ensues: too much cell death can lead to neurodegenerative 



diseases, while too little cell death can lead to autoimmune 
disorders or cancer. Nevertheless, when observing cell death 
under a microscope or plotting a dose-response curve to a 
drug or ligand, one observation recurs frequently, that some 
cells live and others die. How does this cell-to-cell variability 
affect our understanding of cell death? How does it impact 
human health and disease treatments? Although the concepts 
that we aim to cover are likely relevant to all forms of cell death 
and to the treatment of many diseases, we will focus largely on 
apoptosis in cancer, a context in which cell death has been 
heavily studied using a systems biology approach. Triggering 
programmed cell death or circumventing a defect in this 
process have been long-standing goals in cancer therapy and 
the systems biology of apoptosis is so far a pillar of the 
budding field of systems pharmacology. 

How can systems biology help our efforts to understand and 
co-opt cell death in cancer? It may help by extracting 
mechanistic insights in cancer cell behavior and providing a 
broader view of signaling systems, while applying the iterative 
strategy to 'measure, model, and manipulate'. Biologists have 
long applied this strategy, but systems biology makes it more 
explicit by bringing formalized, or mathematical, models to the 
forefront. Here we will discuss how we interpret measure- 
ments of cell death, how cell-to-cell variability impacts it, and 
how we build computational models to understand and predict 
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cell death responses. We also explore the impact of what we 
have learned about the sources of cell-to-cell variability on our 
approaches to manipulating the biochemical state of cancer 
cells to induce cell death. 

Interpreting Measurements of Cell Death in the Presence 
of Cell-to-Cell Variability 

With its emphasis on quantitative results, one of the 
cornerstones of the systems biology of cell death is 
quantitative, or semiquantitative, measurements of various 
steps in the cell death process. Many assays developed to 
characterize programmed cell death since 'apoptosis' was 
coined 1 have been leveraged for systems-level studies. Here 
we consider existing expertly cataloged assays, 2 and the 
information they provide about cell-to-cell variability in cancer 
cell death. 

As Figure 1 illustrates, certain steps in apoptotic cell death 
can be measured by several techniques. These techniques 
vary significantly as to whether or not they can: (1) be 
quantitative, (2) yield single-cell measurements and (3) inform 
us about system dynamics. An illustrative example is to 
consider measurements of caspase activation. One option is 
to detect cleaved caspase substrates from cell lysates 
collected at specific end points by immunoblots. Immunoblots 
can be made semiquantitative - or even quantitative with 
proper calibration curves - and therefore are suitable to a 
systems biology approach (e.g., high-throughput microwes- 
terns). 3 However, because they measure protein abundance 
in homogenates of thousands of cells, by design they do not 
inform on cell-to-cell heterogeneity. With a series of measure- 
ments, one could effectively determine which conditions 
yielded more aggregate caspase activity in the cell population, 
but not how that activity was distributed across single cells 
within this population. 

Another common assay uses fluorogenic caspase sub- 
strates to quantify total caspase activity. 4,5 When used on sets 
of cellular lysates, this assay provides the same kind of 
aggregate, population-level information that immunoblots 
give us. However, if instead one labels intact cells with a 
fluorescent dye that accrues in the nucleus in response to 



caspase activation, the accumulation of caspase activity can 
be measured in single cells, revealing the extent of hetero- 
geneity in the population. 4,5 By assessing a population of 
single cells the experimenter obtains semiquantitative infor- 
mation on caspase activation across this population and can 
determine whether the distribution of caspase activation is 
bimodal (some cells on, some cells off; Figure 2b), or 
unimodal (all cells activating caspases, with some variation 
around the average; Figure 2c). Taking it a step further, 
fluorescent protein-based biosensors designed to detect 
caspase activity enable measurements in individual living 
cells at multiple time points, providing information about the 
dynamics of caspase activation. 6 " 12 Single-cell observations 
of activation dynamics accurately reflect dynamics in a single 
instance of the biochemical system - a cell undergoing a cell 
death decision process - and are therefore often the ultimate 
goal of system biology approaches. 

How much information is missing from population-level 
measurements? That depends on the dynamics and the cell- 
to-cell variability of the process under study and, for the above 
examples, it depends on which caspase is assayed. During 
extrinsic apoptosis, death ligands bind to their receptors and, 
following assembly of a death-inducing signaling complex 
(DISC), activate initiator caspases-8/-10. 13-15 Owing to 
cell-to-cell variations in the abundance of receptors, 
caspase-8, and protein components of the DISC, the timing 
and extent of caspase-8 activation can vary considerably 
between cells exposed to the same death ligand dose 6,11 
Thus, a population-level measurement of caspase-8 activity 
cannot distinguish between a small amount of caspase 
activation in most cells, and a large amount of activation in a 
few cells (Figures 2a and c). In contrast to caspase-8, 
population-level measurements of effector caspase-3 activity 
can effectively report on how many cells have activated the 
protease en route to apoptosis. This is because single-cell 
measurements of caspase-3 activation dynamics have already 
revealed that in extrinsic apoptosis, caspase-3 activation rapidly 
goes from nearly zero to maximal. 6,9 This rapid activation 
results in most cells having either no, or full, caspase-3 
activation at any given time (also observable by flow cytometry; 
Figure 1 and detailed, for example, in Albeck et a/. 6 ). 
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*Noted here as 'Semi-Quantitative' although the method can be made quantitative with the inclusion of appropriate calibration standards. 

Figure 1 Different aspects of apoptosis can be assayed qualitatively or quantitatively at the single-cell level. For each experimental method listed, we note the apoptotic 
cell death processes that it measures as well as whether the measurements can be made in single cells or at the population level, whether they are quantitative or not, and 
whether they enable the tracking of response dynamics. IF, immunofluorescence; IHC, immunohistochemistry; PS, phosphatidylserine 
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Figure 2 Characteristics of heterogeneous and homogeneous signaling responses for different forms of measurements, (a) Schematic representations of population-level 
measurements of signaling and apoptosis over time (left) or for a dose response (right). These are expected to be similar whether the cells respond heterogeneously or 
homogeneously. Dose-response curves can be characterized by the maximal effectiveness (E max ), the half-point of effectiveness (£ 50 ) and the concentration at £ 50 (EC 50 ). 
Recent work (Fallahi-Sichani ef a/. 76 ) has shown that a shallow slope in a dose-response curve (black versus gray dotted line) may indicate heterogeneity in the response, see 
also Box 1. (b and c) Schematic representations of hypothetical results of various single-cell measurements - end point measurements such as flow cytometry (top) and 
dynamic measurements for time-dependent (bottom left) or dose-dependent (bottom right) data series. Although these are all consistent with the hypothetical population-level 
measurements shown in a, they illustrate how population-based measurements can fail to distinguish between heterogeneous all-or-none signals (b) and homogeneous 
graded signals (c) 
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This knowledge of typical single-cell caspase-3 activation 
dynamics should influence our interpretation of population- 
level measurements of extrinsic apoptosis. For example, 
when observing twice as much caspase-3 activity in a sample, 
we can reasonably rule out that caspase activity doubled in 
each cell and instead conclude that there are twice as many 
apoptotic cells. This is a rather context-specific scenario and 
experimenters are still encouraged to use single-cell mea- 
surements, especially when their conclusions depend on 
knowing the distribution of responses within a population. 

Importantly, the single-cell observations of caspase activa- 
tion dynamics, coupled with computational models of caspase 
regulatory networks, have cemented our understanding of 
effector caspase regulation during apoptosis. 6-8,10,16 
Caspase-3 is normally under the tight control of the inhibitory 
protein XIAP (X-linked inhibitor of apoptosis protein). 
However, after mitochondrial outer membrane permeabili- 
zation (MOMP), a key 'point of no return' event in both intrinsic 
and extrinsic apoptosis, 17,18 XIAP is disabled and caspase-9, 
an activator of caspase-3, is activated. This, in effect, is as if 
the cell simultaneously releases the brake (disabling XIAP) 
and presses on the accelerator (activating caspase-9). 
MOMP is therefore an event that promotes the 'snap-action', 
rapid increase in caspase-3 activity characterized in single- 
cell assays. 6,9,12 Computational modeling predicted, and 
subsequent experiments validated, that feedback mecha- 
nisms (e.g., via XIAP cleavage, or caspase-6 activation) are 
not required for the 'snap-action' dynamics. 7 

Implications of Cell-to-Cell Variability for the 
Computational Modeling of Cell Death 

Systems biology models of cell signaling and cellular 
behaviors have been built using several strategies. 19 For 
computational models of both extrinsic (ligand-induced) and 



intrinsic (damage-induced) apoptosis, strategies include 
Boolean or fuzzy logic, 20-22 data-driven 23-30 or Bayesian 
algorithms. 31 Each method is suited to a specific type of 
question, but the modeling approach that can most easily 
relate to the underlying molecular mechanisms is the use of a 
system of ordinary differential equations (ODEs). The system 
of ODEs mathematically encodes the biochemical reactions 
that are understood, or assumed, to take place in a cell death/ 
survival decision. 32,33 Many ODE models of ligand-induced 
cell death have been published, simulating and predicting the 
cell death response to Fas ligand (FasL), TNF-related 
apoptosis-inducing ligand (TRAIL) or tumor necrosis factor 
(TNF) with various degrees of details. 7,10,34-38 

Most ODE models recapitulate a series of biochemical 
reactions that take place over time, and therefore simulate 
response dynamics in a signaling network. It is worth noting 
that one instance of the model represents one instance of the 
biochemical network, and thus represents one cell. Accordingly, 
simulation results should be compared with experimental 
measurements characterizing the dynamic behavior of 
single cells. 

Which cell should a model attempt to recapitulate: a cell 
matching population-averaged measurements or, perhaps, a 
'typical' cell? Although seemingly fastidious, this is a crucial 
decision in the context of models of cell death, principally 
because certain measurable steps of the cell death process, 
like caspase-3 activation, have an 'all-or-none' character, as 
we discussed in the previous section. In general, if there is 
substantial cell-to-cell variability in the cell death process 
being modeled, especially if the cellular responses are 
asynchronous, there may not exist any individual cell within 
the population that responds with a behavior that tracks the 
measured average (Figure 2b). One clear example of the risk 
of relying on population-based measurements when deve- 
loping ODE models of cell death lies with the concept of an 



Box 1 Assessing cell-to-cell heterogeneity using population-based measurements 



The traditional interpretation of a result from a population-based assay is that it defines an expected, or mean, cellular 
behavior. However, the presence of cell-to-cell heterogeneity can still be revealed through careful experimental design and 
thoughtful inspection of population-based data. 



Single-cell heterogeneity from dose-response curves 



One way to determine if cellular responses vary between cells is to closely inspect dose-response curves, which can be 
parameterized with a conventional logistic sigmoid function (equation 1 and Figure Box 1a). 
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Here, / measures the response of cells at dose D (e.g., viability of cells at a given drug concentration), E 0 and E inf are the 
theoretical minimal and maximal effects as defined by the top and bottom asymptotes of the response curve, EC 50 is the 
dose at which the half-maximal effect (E 50 ) is observed and HS is the hill slope coefficient. Although E 0 is measurable, as the 
effect of the absence of drugs, E inf is estimated from the dose-response curve, or from E max , the maximal measured effect. 
For viability measurements, if the theoretical maximum effect (Enf) issuboptimal and therefore does not reach zero or 100% 
(depending on whether the stimulus promotes, or inhibits cell death, respectively), then there must be heterogeneity in the 
cell population. 
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When measuring a signal or an outcome that is not binary or categorical in nature, other more subtle metrics may be more 
telling. Recently, Fallahi-Sichani era/. 76 showed that for the relative viability of cancer cells treated with phosphatidylinositol- 
4,5-bisphosphate 3-kinase (PI3K) pathway inhibitors, a shallow dose-response curve was associated with cell-to-cell 
variability in the degree of signal inhibition. For competitive kinase inhibitors, one would expect the dose-response curve to be 
well fitted with a HS of ~ 1 ; a shallow curve would have HS « 1 (e.g., Figure Box 1). Fallahi-Sichani era/. 76 found that the 
slope of the dose-response curve for some PI3K/Akt/mammalian target of rapamycin pathway inhibitors was particularly 
shallow, with HS « 1 for several cancer cell lines. Single-cell analyses showed that the response to these particular PI3K 
pathway inhibitors had larger coefficients of variation, indicating greater cell-to-cell variability in the population. 



Single-cell heterogeneity from stochastic profiling 

Many experimental techniques are sensitive to the abundance of starting material, limiting their usefulness for single-cell 
analysis. One clever approach to circumvent these limitations, and infer variability within a population of cells, is to make 
several measurements on samples of a few cells (10-20) and quantify the variability between samples. The statistical 
approach, referred to as 'stochastic profiling', reveals heterogeneity, or more specifically bimodality, in the distribution of 
responses across a population of cells. 77 Heterogeneity is detected when the distribution obtained from 12 to 15 replicate 
measurements from small cell populations is broader than expected from simple sampling and measurement error (Figure 
Box 1b). 78 Applying this approach, Bajikar ef al. used microarray gene expression data from 10-cell samples 78 and 
maximum-likelihood inference to reveal a surprisingly large spectrum of single-cell regulatory states in mammary epithelial 
cells in acinar structures. Some of these regulatory states were common (~25% of the cells in each acinar structure), other 
states were rare, (only ~ 1 out of 40 cells); 79 and therefore not previously observed or described. Single-cell heterogeneities 
in gene expression can therefore be deconvolved from population-based experiments by using statistical data models of 
expected measurement distributions. 
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Figure Box 1. Cell-to-cell heterogeneity can be inferred from population-based data, (a) Schematic of drug dose- 
response curves with differing hill slope coefficients {HS). The curves have an increasingly shallow slope as HS becomes 
smaller. Hypothetical data points are plotted for the scenario where HS= 1 (dark blue curve) and the parameters of the 
logistic sigmoid function describing this curve are depicted (£ 0 , E 50 , E max , £j nf and EC 50 ). (b) Hypothetical distributions for a 
homogenously distributed quantitative trait (blue) and for a heterogeneous trait with a bimodal distribution (red) are show at 
the top. At the bottom, a reference, or 'expected', distribution of population-based measurements for the homogeneously 
distributed trait (gray line) is compared with a hypothetical distribution of measurements from small numbers of cells. As the 
proportion of 'high' versus 'low' cells fluctuates from sample to sample, the variability in measurement values is broader than 
the reference distribution if the trait has a bimodal distribution (red histogram), indicating the presence of cell-to-cell 
heterogeneity. 



'£C 50 ' (or half-maximal effective concentration). At the EC S0 
for a certain cell death-inducing drug or ligand, half of the 
treated cells die but, in all likelihood, no cell 'half-dies'. 
Therefore, matching a model solely to a population-average 
measurement can be a foiled exercise and may yield a model 
that does not behave like any cell in the population. 



What if the only accessible quantitative data are from 
population-based assays? As we discuss in Box 1, recent 
work points to ways to analyze population-based data to infer 
whether there is heterogeneity in the cellular responses of 
interest. This could be done by interpreting a dose-response 
curve or by statistical analysis of repeated measurements on 
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small numbers of cells (Box 1 ). Nevertheless, how should we 
then compare model simulation results to population-based 
data? One promising approach to modeling heterogeneous 
cell responses is to run 'populations' of models, to recapitulate 
the behavior of populations of cells. In this approach, sets of 
many simulations of the same model with variations in 
parameters representative of cell-to-cell variability sources 
(the sources of this variability are discussed in the next 
section) are used to approximate the behavior of an ensemble 
of cells. 11,39-44 Averaged simulation results can be compared 
with population-average measurements, or when single-cell 
data are available, sets of simulations can be directly 
compared with sets of single-cell data. Toward this end, the 
recent development of Bayesian and Monte Carlo Markov 
Chain approaches for parameter estimation, where the 
distributions of possible parameter values are derived by 
calculating the best fits to distributions of single-cell measure- 
ments, should prove particularly useful. 45 

In conclusion, when comparing results from simulations and 
experiments, computational modelers should ask: are the 
data semiquantitative, quantitative or qualitative? Do the data 
provide single-cell information? If so, do they provide 
information about single-cell dynamics? Any type of measure- 
ment can be useful, but knowledge of its information content 
will allow the modeler to compare the data with the appropriate 
modeling results. One final point that we have not yet 
addressed is that to accurately predict experimental results, 
we need to understand, with precision, what activity each 
assay measures. Going back to our example of caspase 
activation discussed in the previous section, one must know: 
(1) whether the assay measures the activity of a specific 
caspase or whether it is a measurement of total caspase 
activity (even 'caspase-specific' substrates tend to have 
cross-reactivity with other caspases), 46 and (2) whether the 
assay measures caspase cleavage, caspase activity or 
accumulated caspase cleavage product. For example, it is 
known that cleaved poly (ADP-ribose) polymerase, a product 
of effector caspase activity is relatively stable, and therefore 
can perdure as evidence of past caspase activity. By contrast, 
cleaved caspase-3 itself is relatively unstable, and evidence of 
its activation can disappear over time if monitored solely via 
cleaved caspase-3 abundance. In fact, modeling predicted, 
and experiments validated, that the short half-life of caspase-3 
is important to the dynamics of cell death: if one deletes the E3 
ubiquitin ligase domain of XIAP responsible for tagging 
caspase-3 for degradation, cells can lose the 'snap-action' 
cell death dynamics and may survive with damaged proteins 
and DNA. 8 Just as crucial as comparing population averages 
with simulation averages and comparing single-cell data sets 
to sets of single model instances, accurate pairing of 
measured and modeled entities is key to the development of 
high-quality models. 

Sources of Cell-to-Cell Variability in Cell Death 
Responses 

As we remarked in our introduction, the motivation for single- 
cell approaches to measurements arose from the observation 
that, more often than not, the response to death-inducing 
stimuli is variable. What are the sources of this cell-to-cell 



variability? Which sources are relevant when building systems 
biology models of cancer cell death? 

Perhaps the most obvious source of variation in a 
population of cancer cells is genetic changes. Almost 100 
years ago, the idea of mutations being at the origins of cancer 
first emerged (reviewed in Wunderlich 47 ), and so did the 
observation that cancer cells tend to lose or gain chromo- 
somes. 48 It is well established that chromosomal or mutational 
genomic instability is one of the hallmarks of cancer, an idea 
first proposed nearly 40 years ago (Nowell 49 ; reviewed in 
Negrini etal. 50 ). As a result of this genomic instability, tumors 
are heterogeneous and even cancer cells in culture accumu- 
late mutations, including mutations that lead to anticancer 
drug resistance. For example, treatment with epidermal 
growth factor receptor (EGFR) kinase inhibitors selects for 
cells with mutations in EGFR, 51 ABL (Abelson murine 
leukemia viral oncogene homolog 1) kinase mutations lead 
to resistance of leukemia cells treated with imatinib 52 and 
mutations in the pro-apoptotic protein Bax (B-cell lymphoma 2 
(Bcl-2)-associated X protein) emerge when selecting for cells 
resistant to TRAIL. 53 

When genetic changes directly impinge on the signaling 
events and cellular responses represented in a systems 
biology model, they must be implemented to obtain accurate 
predictions. Different types of mutations require different 
changes to a model: (1 ) if a tumor-suppressor protein species 
is lost, it can be removed from the model; (2) if a mutation 
instead affects binding or enzymatic activity, the appropriate 
reaction rate constants should be modified; (3) if a mutation 
affects protein stability, it can be implemented by changing the 
parameters for synthesis or degradation or, in a less detailed 
model, by changing the concentration of the protein to reflect a 
new equilibrium value. 

But what if one is modeling a cell population heterogeneous 
for a mutation? Then it will become relevant to model a 
population of cells, using a mixture of 'wild-type' and 'mutant' 
models, as we introduced in the previous section. Another 
plausible scenario is that mutations arise as a consequence of 
treatment (as they do for the examples listed above). In that 
scenario, one could model a population of cells and assign to 
each cell a probability for a switch to the mutant phenotype. 
Such an approach has been effectively applied to model the 
evolution of pancreatic tumors and a model prediction, later 
experimentally validated, correctly inferred the effect of 
dosing schedule on the acquisition of resistance in lung 
cancer. 54,55 Efforts such as the Cancer Cell Line Encyclope- 
dia project, aiming to catalog genetic mutations and drug 
sensitivity, will help customize our computational models to 
each experimental model cell line. 56 

Although 'fractional kill', the concept that one round of 
chemotherapy generally does not kill all cells in a tumor, can 
be partly explained by genetic heterogeneity, we now also 
appreciate the importance of non-genetic heterogeneity. A 
common idea for a non-genetic cause of fractional kill is that 
rapidly dividing cancer cells exhibit different drug sensitivity in 
different cell cycle phases. This explanation makes sense for 
cytotoxic drugs that target cell cycle processes - taxanes 
that target microtubules could preferentially act on mitotic 
cells and DNA-damaging agents could be especially toxic for 
cells replicating their DNA or already committed to mitosis. 
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Indeed, rapidly proliferating cells in the body, such as bone 
marrow or hair cells, are specifically targeted by these drugs. 
However, although tumor cells do inappropriately proliferate, 
the fraction of cells in mitosis at any time within a solid tumor 
can be very small, likely too small to explain the amount of cell 
death observed with anti-mitotic drug treatment (reviewed in 
Mitchison 57 and Komlodi-Pasztor et a/. 58 ). Therefore, 
although cell cycle phase impacts the response of cancer 
cells to cytotoxic treatment, other non-genetic sources of 
variability are required to explain the observed response of 
tumors. 

Are other inducers of cancer cell death similarly affected by 
variability in cell cycle phase? For apoptosis-inducing ligands, 
the cell cycle effect hypothesis has been disproved for 
TRAIL, 11 but the question is unresolved for FasL and TNF. 
Early work had reported that TNF lengthened the G2 phase 
and induce death preferentially during or soon after mitosis in 
mouse L929 fibroblasts, 59 while a later study reported that B 
lymphoma cells preferentially die at exit from S phase, and are 
more sensitive to TNF-induced apoptosis if treated during 
G1. 60 The use of single-cell and more quantitative 
approaches should allow us to elucidate whether the cell 
cycle-associated variability in response to TNF is cell type- 
specific or whether there is a generalizable effect. 

When cell cycle effects are observed, or hypothesized, 
accurate models of the cellular responses to apoptosis- 
inducing treatments will have to include them. A rigorous 
approach would be to integrate cell-cycle driving pathways 
into the model, as was done to investigate the interplay 
between cell cycle and DNA-damage response pathways 
regulating cell cycle arrest. 44 This model predicted, and 
experiments confirmed, that knocking out p21 would allow 
endoduplication of DNA in about 50% of y-irradiated cells, 
showing that the model successfully captured the cell cycle 
and DNA-damage response pathway crosstalk 44 A simpler 
approach to integrate cell cycle effects would be to simulate a 
population of models with a data-based distribution of cell 
cycle phases while including, for each cell cycle phase, a 
probability that quantifies the likelihood of pathway activation. 

Cell cycle phase is an example of a potential source of 
variability outside the cell death signaling networks, but what 
about variability within these networks? It is now appreciated 
that the general biochemical 'state' of all cells is somewhat 
plastic, because transcription occurs in stochastic bursts. 61,62 
This ultimately results in variability in the abundance of all 
proteins both across a cell population and in the same cell 
over time (e.g., as observed in Cohen et al. 63 ). Therefore, 
even a deterministic cell death process - which always has 
the same outcome if starting conditions are the same - can 
have a variable outcome because each cell starts with 
different initial protein concentrations. Notably, because the 
state of each cell changes over time, this plasticity actually 
allows cell populations decimated by a death-inducing 
stimulus to repopulate and recapitulate the sensitivity profile 
of the original population after a few days. 11,64 

How can this noise in gene expression be incorporated into 
models of cell death? The most faithful, if cumbersome, 
representation would be to model the synthesis of each 
protein from first principles, with bursts of transcription. 
An effective alternative approach is to model a population of 



cells, as introduced above, by sampling from measured 
distributions of protein concentrations to build a 'population' of 
models. 11,39-42 This approach can be used to mimic any 
measurable, or hypothesized, cell-to-cell variability in protein 
abundance, whether from noise in gene expression or, for 
example, from longer-lasting chromatin-encoded epigenetic 
changes, which have been shown to impact the response of 
genetically equivalent tumor cells to chemotherapeutics 
in vitro 6566 and in xenograft models (Kreso era/. 67 ; reviewed 
in Marusyk and Polyak 68 ). 

In summary, both genetic and non-genetic sources of cell- 
to-cell variability that impact the response to anticancer 
agents have already been identified. Additional sources of 
variability certainly exist in vivo; for example, the microenviron- 
ment differences across different regions of a solid tumor that 
could influence how cancer cells respond to anticancer drugs 
(reviewed in Hanahan and Weinberg 69 ). To design effective 
cancer treatments, we will have to account for these diverse 
forms of variability and incorporate them into our systems 
biology models of cancer treatments. Although this may seem 
like a daunting task, as long as we can measure the source of 
variability we can also model its impact on cellular responses, 
whether the variability arises from outside or inside the 
signaling network of interest. 

Master Regulator versus Multifactorial Regulation - A Road 
to Combinatorial Anticancer Therapies 

Above, we discussed possible sources of cell-to-cell hetero- 
geneity and examples of approaches to account for them in 
systems biology models. But is that complexity necessary? 
What if there was a single master regulator of the cell death 
decision process, a single source of variability? In the 
apoptosis literature, there are countless examples of correlat- 
ing cell fate with the abundance of a key protein. When 
Khaider et al. 70 reviewed factors contributing to TRAIL 
resistance, they found articles citing the importance of 
receptor abundance, receptor trafficking, abundance of 
Bcl-2, Bim (Bcl2-interacting mediator gamma), BID (BH3 
interacting-domain death agonist), PUMA (p53 upregulated 
modulator of apoptosis), Noxa, BAD (Bcl-2 antagonist of cell 
death), Mcl-1 (myeloid leukemia cell sequence 1), Bcl-xl (B-cell 
lymphoma-extra large), FLIP (FLICE (FADD-like IL-1/i-con- 
verting enzyme)-inhibitory protein), STAT5 (signal transducer 
and activator of transcription 5), pro-caspase-8, pro-caspase-3, 
FADD (Fas-associated protein with death domain), p53 and 
even constitutive Akt activation. Which, if any, of these 
cited studies identified the correct 'master regulator' of 
TRAIL-induced cell death versus survival decisions? 
Perhaps tellingly, several implicated two or more regulators. 
Taken together, these studies suggest that the regulation of 
cell death is highly context dependent and multifactorial. 
Although purely experimental approaches often necessarily 
take a reductionist view to identify a key regulator in a specific 
context, systems biology models that formalize our 
accumulated knowledge will enable a multivariate view of 
the system, and ultimately allow reconciliation of many 
experimental findings. 

The multifactorial nature of cell death decisions is a 
powerful rationale for pursuing combinatorial anticancer 
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therapies, where one drug sensitizes cancer cells to the action 
of the other. For example, although TRAIL-based mono- 
therapies have proven unsuccessful, many different co- 
drugging strategies are being investigated (reviewed in 
Hellwig and Rehm 71 ). Essentially, each drug changes the 
biochemical 'state' of each cell over a period of time. This 
'state' is determined by factors such as the relative abundance 
of signaling proteins or cell cycle phase and the values of 
these factors can be used to plot the location of each cell in a 
multidimensional 'state-space'. Co-drugging chemotherapeutic 
strategies aim for one drug to shift cancer cells from a 
state-space region where they are resistant to the second 
drug to one where they are sensitive (Figure 3a). Models and 
experiments can help us map this state-space and predict the 
most promising co-drugging interventions. 

State-space mapping of cellular responses to death- 
inducing stimuli will necessitate predictive, well-validated 
models. For ligand-induced apoptosis, high-quality models 
exist and mapping results are emerging. For example, 
Howells et al. 72 created a model to identify which cellular 
parameters determine whether abundance and phosphoryla- 
tion of BAD, a pro-apoptotic Bcl-2 family protein, influence 
MOMP and ultimately cell fate. Although not yet experimen- 
tally validated, this model helps predict scenarios where BAD- 
mimicking drugs would increase sensitivity to death-inducing 
ligands. In this study, a cellular response map was derived 



from bifurcation analysis, a mathematical method that maps 
where steady-state cellular responses transition from one to 
another (e.g., to map the blue plane separating 'resistant' and 
'sensitive' states in Figure 3a). However, not all apoptosis- 
inducing signals evolve to a steady state. When transient 
events dominate the cell death decision, bifurcation analysis 
becomes inadequate. For the response of cells to CD95 
(cluster of differentiation 95, also known as Fas receptor) 
activation, Neumann et al. 36 aptly used forward model 
simulations, predicting system behavior within a multidimen- 
sional grid of initial conditions, to map regions of the cellular 
FLIP (c-FLIP) versus pro-caspase-8 space that lead to strong 
activation of pro-survival (NF-/vB) nuclear factor kappa-light- 
chain-enhancer of activated B cells signaling versus pro- 
apoptotic caspase-3. This map predicted how to manipulate 
c-FLI P and/or caspase-8 abundance or activity to promote cell 
death, and these predictions were confirmed experimentally 
using genetic manipulations to vary protein abundance and 
small molecule inhibitors of caspases to manipulate caspase 
activity. 36 

A third method, the calculation of Lyapunov exponents, can 
similarly be applied to highlight where system (or cellular) 
behavior diverges in state-space (see also overview in Box 2 
of Aldridge et a/. 8 ). Lyapunov exponents were applied to a 
seven-dimensional state-space of a computational model of 
TRAIL-induced cell death and predicted that the cellular 





Figure 3 Multivariate and dynamic analyses of cell state to predict synergistic effects of drug regimens on cellular responses, (a) Hypothetical three-parameter 'state- 
space' maps showing the position of a population of cells (blue), as well as the plane separating regions of state-space leading to sensitivity (back) or resistance (front) to 
treatment. Here 'parameters' could quantify the abundance or activity of certain proteins, for example. In the scenario illustrated here, the values of parameters 1 and 2 are 
found to govern the sensitivity of the cells to a given drug. The middle graph shows the results of a moderately successful intervention, which results in many cells becoming 
sensitive to the drug. The graph on the right shows a successful intervention where the entire population becomes sensitive. Note that interventions can also lead to changes in 
other parameters, such as parameter 3, which do not affect the sensitivity of the population, (b) Graphs showing hypothetical cell death response surfaces for treatment with 
various concentrations of drugs 1 and 2. The time interval between administration of drug 1 and drug 2 increases to the right, potentially leading to much stronger or weaker 
effects than the simple additive effect achieved by the simultaneous administration seen in the leftmost panel. The vertical axis (death (%)) and the color scale (from blue (0%), 
to cyan, green, yellow, orange and red (100%)) both indicate the fraction of cells dying in the population 
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concentration of several proteins, namely ligand, receptor, 
initiator caspases-8 (and -1 0), XI AP and caspase-3 influenced 
whether or not cells required MOMP to commit to extrinsic 
apoptosis. 8 Experiments confirmed that, indeed, XIAP to 
caspase-3 ratio and ligand concentration were major deter- 
minants of the requirement for MOMP. Interestingly, another 
prediction from this analysis was that while in certain 
biochemical contexts (or regions of state-space), cell-to-cell 
variability in XIAP and caspase-3 concentration should cause 
phenotypic cell-to-cell variability, in other contexts all cells 
were predicted to have the same phenotype. This was 
validated experimentally by showing that whereas the fraction 
of MOMP-dependent cells in the T47D breast carcinoma line 
was ligand concentration dependent, other cell lines, such as 
the SKW6.4 B-cell lymphoma line, had no phenotypic 
heterogeneity and were even insensitive to overexpression 
of XIAP. 8 This context dependence of the impact of variability 
in protein abundance on cellular responses likely arises in part 
from the underlying network topology and can be explored 
using state-space maps (see also Gaudet et al. 39 ). 

Using these mapping approaches and others, we can start 
to understand how to maximize the impact of co-drugging 
strategies by re-positioning cancer cells into region of the 
'state-space' where the entire population is drug sensitive 
(Figure 3a). As discussed above, in concept, drug treatments 
shift the biochemical 'state' of a cell. This shift is dynamic and 
thus the position of cells will depend on time since drug 
addition. Lee et al. 23 expertly demonstrated this principle, a 
process they refer to as dynamic network rewiring. They 
optimized relative timing and sequence of drug addition 
(Figure 3b); over time, the first drug 'rewires' the signaling 
network, by moving cells within state-space, maximizing 
responses to subsequent drugs. They found that treatment 
with the EGFR inhibitor erlotinib 4h before addition of 
doxorubicin, a conventional DNA-damaging chemotherapy, 
markedly enhanced killing of triple-negative breast cancer 
cells compared with either drug alone, both drugs added 
simultaneously, or doxorubicin preceding erlotinib. 23 A data- 
driven model was built from systematic time-dependent 
measurements within the responsive signaling networks 
following EGFR inhibition by erlotinib, and the model identified 
a key predictive dynamic biomarker for cell sensitization. Late 
increase in cleaved caspase-8 was specifically associated 
with synergistic effects of sequential erlotinib-doxorubicin 
treatment, likely indicating reactivation of otherwise dysregu- 
lated extrinsic apoptosis. 23 Although no simple systematic 
way exists to probe whether synergy could arise from specific 
sequential or co-drugging regimens, it is clearly a promising 
frontier for the systems biology of anticancer treatments. 

The shape of cleaved caspase-8 accumulation over time is 
a biomarker for the synergistic response to sequential erlotinib 
and doxorubicin treatment, and there are other examples 
where signaling dynamics encode information for cellular 
responses. One striking example involves the p53 tumor 
suppressor: in single cells, y-irradiation induces oscillations in 
p53 nuclear abundance and cell-cycle arrest, while ultraviolet 
treatment induces sustained p53 nuclear localization and 
apoptosis (reviewed in Purvis and Lahav 73 ). Could pharma- 
cological interventions that manipulate signaling dynamics be 
used to induce specific cell fates? Purvis era/. 74 showed that 



this could be done for p53: as predicted by their model, timed 
treatments with nutlin-3, an inhibitor of MDM2 (mouse double 
minute 2 homolog)-directed degradation of p53, forced 
sustained elevated nuclear abundance of p53 after y- 
irradiation, resulting in cell death. Such manipulations could 
also apply to other well-characterized systems, as enticingly 
proposed by Behar et al. 75 They used computational models 
to virtually screen for conditions yielding specific dynamics, 
then applied their findings to NF-/cB-driven transcription. They 
found that, as predicted by their analysis, different inhibitor 
treatments yielded stimulus-specific effects on early versus 
late NF-kB target gene expression. 75 By combining accurate 
measurements of cell 'states' and single-cell response 
dynamics, we can learn how signaling dynamics encode cell 
fate, and use models to predict how to manipulate these 
dynamics to obtain the desired cellular outcomes. 



Conclusion 

There is an oft-quoted phrase from George EP Box, a British 
statistician, stating that 'all models are wrong, but some are 
useful'. As models are by necessity an approximation of 
reality, all models are 'wrong'. However, when carefully 
developed and solidly anchored in high-quality data, models 
can be predictive, leading to testable and falsifiable hypoth- 
eses and thereby allowing us to learn about the system under 
study. Here we have discussed the frequently observed 
phenotypic heterogeneity in the response of human cells to 
various death-inducing stimuli as well as the sources of this 
variability and implications for measuring and modeling cell 
death. We are entering an era of 'biology in the second 
moment' where the focus is not the average, or dominant, 
behavior exhibited by a population of cells, but rather the focus 
is on the variance, and we anticipate that this new focus will 
continue to contribute to our understanding of the regulation of 
cell death in health and disease. 
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